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I.  INTRODUCTION 


The  elliptical  coverage  function  P(R,  <rx,  <7y,  h,  k)  represents  the  probability  of  a  shot  or  event 
occurring  within  a  circle  C  in  the  xy-plane  with  radius  R  and  center  (h,k),  i.e., 

C:  (x-h)2  +  (y-k)2  =  R2  . 

The  shot  falls  under  an  uncorrelated  bivaria'e  normal  distribution  with  mean  (0,0)  and  standard 
deviations  ?x,  <?y  in  the  x  and  y  directions,  respectively.  Thus 

P(R,  ax,  Wy,  h,  k)  =  --  f  |  E[Xj/(^  ax)]  E{yt/(^2  ay)}  dx:  dyx  (1) 

X  y  *  Area  of  C 

P(R,  Wx,  Wy,  h,  k)  =  P(R,  ^y,  ffx,  k,  h)  =  P(P,  ffx,  ay,  |h|,  |k|),  (2) 


E(z)  =  exp(--z2). 


Letting  Xj  =  ^2  ax  x,  and  yt  =  ^2  cry  y,  one  obtains 

( h  +  R)/(>12  ?x)  [k  +  >J  R  2  -  (h  - >(2  ax  x) 2  j  / (>f2  Wy) 

P  =  l|  Efx)  |  E(y)dydx.  (■ 

(h-R)/(^ax)  [k_>|R2_^_^-xX)2j  j (^Sy) 

It  is  easily  seen  from  (4)  that  P  is  a  function  of  four  independent  variables.  Therefore,  taking 
(2)  into  account,  we  assume  hereafter  that  ?x  >  ay  >  0,  and  that  all  variables  are  normalized  by  <xx 
with  an  additional  factor  of  l/-\f2  for  R,  h,  k.  Hence 

R~W^'  h  =  ^jA^’  k  =  V^’  ay  =  ^(^1)’  <Tx=1-  (; 

Splitting  the  x-integration  interval  in  (4)  into  [h  -  R,  h]  and  [h,  h  +  R]  and  using  the 
transformation  x  =  h  —  Ru  in  the  first  integral  and  x  =  h  -|-  Ru  in  the  second  integral  gives 

P  =  2^=Jp  [E(h-Ru)  +  E(h  +  Ru)jaerf(ky,  Ry-^1  -u2)  du  ,  (i 


ky  k  j  (T y ,  Ry  R/cZy 


and,  for  all  real  p  and  q, 
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aerf  (p,  q)  =  erf  (p  +  q)  -  erf  (p  -  q)  (8) 

erfpH?w|o  dz' 

In  (6),  the  square  root  singularity  at  u  =  1  is  eliminated  by  letting 

t  =  ’Jl  -  u  . 

Hence, 

P=^|b{E[h-R(l-t2)]  +  E[h  +  R(l-t2)]}aerf(ky,  Ryt^-t2)  tdt  (9) 

where  a  =  0  and  b  —  1  . 

Relation  (9)  appears  in  [2,  3,  4].  In  those  papers  the  objective  was  to  give  a  numerical  pro¬ 
cedure  to  compute  P  to  6  decimal  digit  absolute  accuracy  with  the  practical  limitations: 

1/15  <  <  15,  0  <  h  <  600,  0  < >12  k  <  600,  10"6  <  P  <  1  -  10 "6  .  (10) 

In  this  report  an  algorithm  is  given  for  computing  P  which  is  based  on  relative  rather  than  absolute 
error.  A  new  subroutine,  written  in  Fortran,  using  this  algorithm  is  contained  in  the  Naval  Surface 
Warfare  Center  (NAVSWC)  Library  of  Mathematics  Subroutines  (NSWCLIB),  [8].  The  new  subroutine 
is  called  PRILL  and  it  replaces  the  previous  version  which  was  also  called  PRILL.  P  is  now  computed 
to  at  least  6  significant  digits  provided  that 

10"  20  <  P  <  1  -  10  "10,  and  max[h,  ky,  1/0J5  <ry)]  <  5  -J2/EPS  , 

where  EPS  denotes  the  smallest,  machine  dependent,  positive  number  such  that 

1.0  + EPS  >1.0  .  (11) 

For  the  IBM  PC  double  precision  arithmetic  EPS  =  2 "  52  ~  2.222  10  " 16  ;  for  the  CDC  6000  -  7000 
series  single  precision  arithmetic  EPS  =  2"  47  ~  7.105  10 -15.  If  <?y,  h,  k  and  P  satisfy  (10),  then  P  is 
accurate  to  at  least  8  significant  digits. 

It  is  worth  mentioning  that  there  is  much  to  be  gained  from  the  increased  robustness  and 
accuracy  made  available  by  the  new  subroutine,  PRILL.  Indeed,  in  the  future  it  will  be  possible  to 
design  an  inverse  procedure  by  which  R  is  found  as  a  function  of  P,  Oy,  h,  k.  Also,  since  P  is  used  to 
compute  its  3-dimensional  analog  [5],  increased  robustness  and  relative  accuracy  will  also  be  achieved 
for  that  computation. 
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n.  DISCUSSION  OF  THE  ALGORITHM 


In  order  to  compute  P  with  relative  accuracy  using  (9),  or  its  equivalent  with  the  order  of 
integration  reversed,  three  main  sources  for  loss  of  digits  are  recognized: 

(a)  If  try  <C  1,  inherent  error  can  play  a  major  role  when  R  is  large  and  is  very  close  to  k  oi 
\fh2  +  k2  .  When  this  inherent  error  dominates  the  generated  error,  PKILL  does  the  best  it 
can  with  no  claim  on  relative  accuracy. 

(b)  A  loss  of  digits  was  possible  in  the  past  in  the  computation  of  aerf(x,  y)  when  y<x.  A 
routine  has  been  written,  called  AERF  and  contained  in  NSWCLIB,  which  gives  aerf  to  12 
significant  digits.  It  is  described  in  Appendix  A. 

(c)  Truncation  error  can  occur  due  to  using  32  order  Gaussian  quadrature  to  numerically 
evaluate  P. 

Extensive  testing  was  carried  out  to  insure  the  accuracy  claimed.  This  testing,  using  a  Compaq 
Deskpro  386/20  PC,  compared  PKILL  double  precision  results,  with  the  results  from  a  double  precision 
adaptive  quadrature  subroutine,  DQAGS,  which  is  contained  in  NSWCLIB.  Also,  some  additional 
checks  were  carried  out  by  Alfred  Morris,  Jr.  (NAV^V"  }  ..  Tg  an  integrator  he  developed  for  this 
purpose. 

PKILL  initially  processes  a  few  tests  to  check  if  P  can  be  set  to  1  or  0.  P  is  set  to  0  when  these 
tests  indicate  that  P  <  10  “  50  ;  if  P  is  estimated  to  be  greater  than  1  —  10  ~ 10  then  P  is  set  to  1.  These 
tests  are  described  in  [4]. 

If  <Ty  is  sufficiently  small  with  R  >  k,  then  P  can  often  be  evaluated  without  resorting  to 
numerical  quadrature.  We  reverse  the  integration  in  (4)  and  use  (5)  to  get 


wheie 


ky  +  Ry 

P-if  E(y)  J(y,  <ry)  dy  , 
ky-Ry 


(12) 


For  convenience,  J(y,<7y)  will  also  be  denoted  by  J. 


If  R  >  k,  then  for  sufficiently  small  <Ty ,  the  integration  limits  on  y  in  (12)  can  be  set  to  -oo 
and  oo.  We  proceed  by  obtaining  the  Maclaurin  expansi  m  of  J  in  <7y.  The  expansion  to  second  order  is 


where 

J  ~  J0  +  J0'  <7y  +  J0"  <7y2/2  ,  J0  =  J(y,  0),  J0'  =  ^(y,  0)  , 

(13) 

J0  =  -j-  aerf(  h,  >|R2  -  k2  ) 

(14) 
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V  =-r— r[E(h"  ^R2~k2)  +E(h+  >lR2-k2)] 

NR  -k 


r?{(-R^+#rp-2lJ)E(h- ,li?::F) 


(15) 


it  _ _ y* 


■Jr? 


^s6?+^5+2iJ)E0,+  '1Fri5)}' 


(16) 


f  00  f  oo 

Now  E(y)  J0'  dy  =  0  since  y  E(y)  dy  =  0.  Therefore  for  small  <ry,  P  can  be  approximated 

J  -oo  J  -oo  J 


by 


provided  that 


P  ~  ^E(y)  J0  dy  -  ^  J0  =  i  aerf(  h,  >jR2-k2  ) 

<ry2  £  E(y)  |J0"|  dy  <  e  aerf(  h,  -jR2-k2). 

J  J  -oo 


(17) 

(18) 


In  PKILL,  (17)  is  used  when 

W  {(  ^  k2  +  2k2l  h  ->|R2  -  k2  1 )  E(h  -  >|R2  - k2  )  |  <  c  aerf(  h,  ^R2-k2).  (19) 

The  left  hand  side  c.  (19)  is  an  upper  bound  for  the  left  hand  side  of  (18)  since 

H  y2E(y)dy=f . 

J  -00  L 

The  quantity  e  is  given  by 

e  =  min  (6.71  ^JEPS,  10  ~  5 ).  (20) 

For  the  case  of  sufficiently  small  <ry  and  R  =  k,  a  Similar  analysis  is  used.  Now 

J  =  aerial-,  ~z  ^1  -  cz2)  (21) 

where  _ 

r  =  ^2Ry  ,  c  =  y/(2R),  z  —  ^py  . 

Expanding  J  to  third  order  in  z  gives 

J  ~  J0  4  J0'  z  +  J0"  z2/2  4  V"  z3/  6-  (22) 

But  J0  =  0  and  J0"  =  0,  so  that 

J~Jo'z4-Jo'"z3/0-  (23) 

Then  from 
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|I  =  2  C  [E( h  +  r  z  >|l  —  c  z 2  )  +  E( h  -  r  2  -\jl  —  c  z2  )]  , 


we  have 


and  from  c^J/dz3  we  obtain 


V  =  2  r  E(h)  , 


J0"76  =  2  r  [  -  c/2  +  r2  (2  h2  - 1)/3]  E(h). 


Hence 


provided  that 


where 


P~2z>[2RE(h)l|^0E(y)^y  dy 

P~  |j2R^E(h)J^E(y)^-dy  =  ^E(>>)  r(|)  J2R^ 

2<ry|  |R(2h2_1)  +  4t  |  j0  E(y) y \|ydy / r(|)  <e, 
E(y)  y  JT  dy  =  r(| )  /2,  r(|)/r(|)  ~  .73966878. 


If  the  approximations  for  P  do  not  apply  then  (9),  or  the  corresponding  relation  with  the  order 
of  integration  reversed,  is  evaluated  by  using  a  32  order  Gaussian  quadrature.  The  Gaussian  abscissas 
and  weights  are  denoted  by  V(|j|)  and  W(|j|),  respectively,  jj|  <  16,  j  ^  0,  [1,  p.917].  Thus  from  (9), 
with  b  and  a  denoting  the  upper  and  lower  limits  of  integration,  P  is  approximated  by 

33 

w  £  w(tj|)GjV  j=n-17’  <3°) 

n  =  1 
n?£17 

where 

t.  =  signO)^V(|j|)  +  ^,  (31) 

z(  tj)  =  Ry  tj  ^  2  —  tj2  , 

Gj  =  {E[h  -  R (1  -  tj2 )]  +  E[h  +  R (1  -  tj2 )]}  aerf [ky,  z  ( tj  )]  . 

The  length  of  the  integration  interval,  b-a,  can  often  be  reduced  because  of  the  spike-like 
character  of  the  integrand  in  (9)  for  large  h  and/or  k.  In  [3,  p.378],  a  parameter  a  was  introduced 
which,  for  a  given  c,  led  to  an  effective  integration  interval  [a,b]C[0,l].  That  determination  of  a 
was  based  on  absolute  accuracy  requirements.  Here  this  parameter,  which  we  now  call  A,  depends  on 
the  value  of  P  because  of  the  more  stringent  relative  accuracy  demands.  It  ranges  from  6.5  to  16. 

Rough  order  of  magnitude  estimates  of  P  are  used  to  determine  A.  These  estimates  are  given 
by  T9,  YY,  and  ZZ  below. 
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f  R  E(h)  aerf(ky,  R^)  Oy  <  .05, 
l  1  otherwise 

YY  =  |  aerf(h,  R)  aerf(ky,  Ry)  >  P 
ZZ  =  7(c,  x)  (Incomplete  gamma  function, 
c  =  [l+<ry2  +  2(h2  +  k2)]2/T 

x  =  2  R2  [l  +  <Ty2  +  2  (h2  +  k2)]/T 

T  =  2  [l  +  cry4  +  4  (h2  +  k2  <ry2)]  . 


ky>Ry 


[6]) 


(32) 


The  expression  ZZ  was  given  by 
by  numerical  testing,  is  denoted 


F.  Grubbs,  [7].  The  estimate  we  choose  for  P,  which  was  determined 
by  Zj  and  is  given  by 

[YY  if  YY  <  5-10  - 15 

[min  ( YY,  ZZ )  otherwise  . 


={: 


The  procedure  for  selecting  A,  using  Zj  and  T9,  is  given  in  the  Fortran  listing  of  PKILL  in  Appendix  B 
(see  page  B-5 ). 


It  was  shown  in  [2,  3]  that  for  a  given  value  of  A,  with  /?  =  A/^2,  the  interval  of  integration 
may  be  reduced  according  to  the  following  expressions: 


b  =  < 

. 


gi-(h-/J)/R 

1 

0 


if  0  <  (h  —  /?)/R  <  1 
if  (h-/?)  <  0 
if  (h-/?)/R>  1 


^-(h  +  ^/R  if 
0  if 


(h  +  0)/R<l 
(h  +  /?)/R  >  1  . 


(33) 


(24) 


The  interval  of  integration  may  possibly  be  reduced  further  by  the  following  considerations.  Let 

f(t)  =  aerf[ky,  z(t)j,  z(t)  =  Ryt^-t2  (see  (9))  . 

The  function  ky  -  z(t)  decreases  from  ky  to  ky  -  Ry.  We  already  know  that  if  ky  -  Ry  >  /?,  then  f  ~  0 
since  6.5/>|2  <  /?  <  16/^2.  Hence  there  may  exist  a  value  of  t,  say  t  ,  such  that  ky  -  z( t)  =  /?  in  which 
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case  whe  integration  in  (30)  need  only  be  carried  out  from  a  =  max[t ,  a  (from  (34)]  to  b.  The  quantity 
t  is  found  explicitly  from  ky  -  z(  t )  =  0  and  is  given  by 

t  =  a/[l+^l  -a2] l!2  ,  provided  0  <  a  <  1,  where  a  =  (ky  -  /?)  /  Ry.  (35) 

Now  let 

E  =  (b  -  a)/2,  D  =  (b  +  a)/2,  F  =  l-D, 

where  E  and  D  appear  in  (30)  and  (31)  and  F  is  used  in  (30)  to  compute  1  - ty2  in  Gy.  In  order  to  sum¬ 
marize  the  results  for  E,  D,  and  F  let 

H  =  (h  —  /?)/R,  Z  =  (h  +  /?)/R,  U=n|1-«2  .  (36) 

We  note  that  H  <  1.0,  otherwise  P  ~  0.0  as  noted  in  (33). 


Then  if  H  <  0,  Z>1, 


If  H  <  0,  Z  <  1 , 


D  =  (l  +  a/>f(ITu)/2,  E-  F  =  U/(4D),  ifU<Z, 
E  =  D  =  F  =  1/2,  otherwise. 

D  =  (l  +  a/4(rTU)/2,  E=F  =  U/(4D),  if  U  <  Z, 
D  =  (l+Vn^Z)/2,  E  =F  =  Z/(4D),  otherwise. 


If  H  >  0,  Z  <  1  , 

D  =  [>rrrH  +  >[rru]/2,  E  =  (U-H)/(4D),  2F  =  r^1--  +  .-..^— ,,  if  H  <  U  <  Z, 

D  =  HT^ +  >fT^Z]/2,  E  =  13/(2  R  D),  2  F  =  — B  _  =  +  — ,  otherwise. 

If  H  >  0,  Z>1, 

D  =  [<T=E  +  jT=V]/2,  E  =  (U  -  H)/(4  D),  2  F  =  +  T+^V  '  ifH<U<Z’ 

E  =  D  =  2>J1'-H  ,  F  =  ^  (1  +^—J===  ) ,  otherwise. 


Once  E,  D,  and  F  are  computed,  the  order  of  integration  is  reversed,  by  interchanging  h,  k  e  nd 
irx,  Wy,  and  E,  D,  and  F  are  computed  again.  Comparisons  of  the  two  sets  of  results  for  E,  D,  and  F 
are  then  used  to  determine  the  order  of  integration  to  be  used  in  evaluating  P.  The  details  can  be  seen 
by  looking  at  the  PKILL  code  in  Appendix  B  (see  page  B-8). 
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m.  NUMERICAL  RESULTS 

The  accuracy  of  the  results  from  PKILL  was  established  by  extensive  computer  testing.  An 
inverse  procedure,  using  PKILL,  gave  values  of  R  for  given  values  of  P0,  <rx,  cry,  h,  k .  Then  an 
adaptive  integration  subroutine,  called  DQAGS,  which  is  contained  in  NSWCLIB  [8],  was  used  to 
compute  values  of  P,  denoted  by  Pj,  for  given  values  of  P0,  R,  7fx,  Wy,  h,  k.  The  results  were 
compared  with  those  from  PKILL,  and  the  relative  error  (RERR)  was  found  from  |P^  -  Pj^I/Pq, 
where  Pj^  was  obtained  from  PKILL.  A  few  numerical  results  are  tabulated  below. 

P0  =  1  •  10  ~~ 10  gx  =  1.0 


Case  No. 

R 

h 

k 

CTy 

RERR 

1. 

1.00000000003092E-5 

0.00 

0.00 

.500 

0.0E-00 

2. 

4.47213595556272E-6 

0.00 

0.00 

.100 

9.0E-16 

3. 

4.41421356414086E-6 

0.00 

0.00 

.010 

2.6E-16 

4. 

2.71828182707694E-5 

0.00 

1.00 

.500 

2.6E-16 

5. 

4.04044149 100285E- 1 

0.00 

1.00 

.100 

2.6E-11 

6. 

9.41563620089423E-1 

0.00 

1.00 

.010 

1.7E-14 

7. 

1.28402541674078E-5 

1.00 

0.00 

.500 

1.2E-14 

8. 

5.74233623428373E-6 

1.00 

0.00 

.100 

7.0E-15 

9. 

1.81588616245426E-6 

1.00 

0.00 

.010 

5.3E-15 

10. 

1.93372255656656E0 

1.00 

5.00 

.500 

4.5E-15 

11. 

4.39270699558773E0 

1.00 

5.00 

.100 

1.1E-09 

12. 

4.94102329874250E0 

1.00 

5.00 

.010 

2.4E-14 

13. 

1.40747750223962E-2 

5.00 

1.00 

.500 

3.9E-16 

14. 

6.432548 1 1832100E- 1 

5.00 

1.00 

.100 

1.3E-16 

15. 

9.67161512902738E-1 

5.00 

1.00 

.010 

5.2E-16 

16. 

7.50347341230299E0 

5.00 

10.0 

.500 

2.6E-10 

17. 

9.56989445416763E0 

5.00 

10.0 

.100 

2.0E-12 

18. 

9.96348038065810EO 

5.00 

10.0 

.010 

6.5E-16 

19. 

5.6 199 1593272086E0 

10.0 

5.00 

.500 

2.8E-12 

20. 

6.15159059893341E0 

10.0 

5.00 

.100 

2.7E-09 

21. 

6. 18350703594593E0 

10.0 

5.00 

.010 

1.1E-08 

22. 

5.06546585 192287E2 

100. 

500. 

.500 

5.8E-09 

23. 

5.08536986866509E2 

100. 

500. 

.100 

1.4E-09 

24. 

5.08690971362607E2 

100. 

500. 

.010 

4.7E-08 
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P0  =  1-1Q~10  <rx=1.0 


Case  No. 

R 

h 

k 

<ry 

RERR 

25. 

5.03634 144947275E2 

500. 

100. 

.500 

8.4E-08 

26. 

5.03664419834487E2 

500. 

100. 

.100 

6.5E-08 

27. 

5.0366567  5807472E2 

500. 

100. 

.010 

8.9E-08 

28. 

9.99997319325154E5 

1000. 

1E6 

.500 

2.0E-08 

29. 

9.99999863834799E5 

1000. 

1E6 

.100 

2.2E-08 

30. 

1.00000043606992E6 

1000. 

1E6 

.010 

2.1E-07 

31. 

9.99994138661486E5 

1E6 

1000. 

.500 

2.1E-08 

32. 

9.99994138662153E5 

1E6 

1000. 

.100 

1.7E-07 

33. 

9.99994138662153E5 

1E6 

1000. 

.010 

1.5E-09 

P0  =  .0050 

=  1.0 

Case  No. 

R 

h 

k 

ay 

RERR 

34. 

7.082 14973500709E-2 

0.00 

0.00 

.500 

8.7E-16 

35. 

3.rB243343370606E-2 

0.00 

0.00 

.100 

1.4E-15 

36. 

1.06881353713506E-2 

0.00 

0.00 

.010 

6.9E-16 

37. 

1.87773801239816E-1 

0.00 

1.00 

.500 

5.2E-16 

38. 

8.10127444936150E-1 

0.00 

1.00 

.100 

8.0E-13 

39. 

9.85554779402774E- 1 

0.00 

1.00 

.010 

6.6E-11 

40. 

9.09820386010898E-2 

1.00 

0.00 

.500 

1.7E-15 

41. 

4.10294072289614E-2 

1.00 

0.00 

.100 

5.2E-16 

42. 

1.43927001021222E-2 

1.00 

0.00 

.010 

1.0E-15 

43. 

3.84954308 155358E0 

1.00 

5.00 

.500 

1.0E-13 

44. 

4.79375448729158E0 

1.00 

5.00 

.100 

1.5E-12 

45. 

4.98408394392712EO 

1.00 

5.00 

.010 

9.1E-11 

46. 

2.62868451641037E0 

5.00 

1.00 

.500 

1.7E-16 

47. 

2.62253474480775E0 

5.00 

1.00 

.100 

2.4E-12 

48. 

2.62232994242921  E0 

5.00 

1.00 

.010 

2.6E-12 

49. 

9.64438962680617EO 

5.00 

10.0 

.500 

5.8E-13 

50. 

1.02483865 120650E1 

5.00 

10.0 

.100 

2.4E-12 

51.. 

1.02892071201973E1 

5.00 

10.0 

.010 

3.2E-12 
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P0  =  .0050  <7X  =  1.0 


Case  No. 

R 

h 

k 

RERR 

52. 

8.85881461639006E0 

10.0 

5.00 

.500 

1.3E-11 

53 

8.94686556206516E0 

10.0 

5.00 

.100 

4.3E-12 

54. 

8.95084289534896E0 

10.0 

5.00 

.010 

0  6E-12 

55. 

5.08543123356466E2 

100 

500 

.500 

•1.6E-10 

56. 

5.09342260944039E2 

100 

500 

.100 

1.9E-10 

57. 

5.09402405988742E2 

100 

500 

.010 

3.6E-12 

58- 

5.0  /363938745357E2 

500 

100 

.500 

7.8E-12 

59. 

5.07375894420008E2 

500 

100 

.100 

2.8E-12 

60. 

5.073763893982 12E2 

500 

100 

.010 

2.3E-12 

61. 

9.99999212083791E5 

1000 

1E6 

.500 

2.2E-11 

62. 

1.00000024240470E6 

1000 

1E6 

.100 

5.4E-10 

63. 

1.0000004741 1365E6 

1000 

1E6 

.010 

1.6E-07 

64 

9.99997924171662E5 

1E6 

1000 

.500 

8.2E-11 

65. 

9.99997924171859E5 

1E6 

1000 

.100 

2.2E-08 

66. 

9.99997924171859E5 

1E6 

1000 

.010 

2.7E-11 

P0  =  .500 

°-x 

=  1.0 

Case  No. 

R 

h 

k 

RERR 

67. 

8.70417428244252E-1 

0.00 

0.00 

.500 

0.0E-00 

68. 

6.81985088272086E-1 

0.00 

0.00 

.100 

6.7E-16 

69. 

6.74563888092905E-1 

0.00 

0.00 

.010 

2.5E-12 

70. 

1.35499509964184E0 

0.00 

1.00 

.500 

1.3E-15 

71. 

1.22603637241642E0 

0.00 

1.00 

.100 

4.3E-15 

72. 

1.20638167322879E0 

0.00 

1.00 

.010 

2.5E-12 

73. 

1.18221140061131E0 

1.00 

0.00 

.500 

6.7E-16 

74. 

1.05532095681235E0 

1.00 

0.00 

.100 

1.9E-12 

75. 

1.05059188896165E0 

1.00 

0.00 

.010 

2.6E-12 

10 
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P0  =  .500 

<7X 

=  1.0 

Case  No. 

R 

h 

k 

(Ty 

RERR 

76. 

5. 18608676231 1 12E0 

1.00 

5.00 

.500 

8.9E-16 

77. 

5.13530939278725E0 

1.00 

5.00 

.100 

2.9E-15 

78. 

5. 10946764748624E0 

1.00 

5.00 

.010 

2.2E-12 

79. 

5.12431831343601E0 

5.00 

1.00 

.500 

1.0E-12 

80. 

5.10003905042201EO 

5.00 

1.00 

.100 

7.7E-11 

81. 

5.09902971 170382E0 

5.00 

1.00 

.010 

8.0E-11 

82. 

1.12084500495804E1 

5.00 

10.0 

.500 

4.2E-15 

83. 

1. 1 1824977123348E1 

5.00 

10.0 

.100 

2.4E-12 

84. 

1.11803622401548E1 

5.00 

10.0 

.010 

2.6E-12 

85. 

1.1 1935071779149E1 

10.0 

5.00 

.500 

6.9E-11 

86. 

1.1 1808975580793E1 

10.0 

5.00 

.100 

7.9E-11 

87. 

1.1 1803454776247E1 

10.0 

5.00 

.010 

8.0E-11 

88. 

5.09902830500163E2 

100. 

500. 

.500 

8.0E-11 

89. 

5.09902155326263E2 

100. 

500. 

.100 

2.4E-12 

90. 

5.09901953902433E2 

100. 

500. 

.010 

3.4E-12 

91. 

5.09902203786139E2 

500. 

100. 

.500 

8.0E-11 

92. 

5.09901961553339E2 

500. 

100. 

.100 

8.0E-11 

93. 

5.09901951359279E2 

500. 

100. 

.010 

8.3E-08 

94. 

1.00000050000037E6 

1000 

1E6 

.500 

1.8E-10 

95. 

1.00000050000037E6 

1000 

1E6 

.10 

4.1E-10 

96. 

1.00000050000037E6 

1000 

1E6 

.010 

1.3E-10 

97. 

1.00000049999987E6 

1E6 

1000 

.500 

1.0E-07 

98. 

1.00000049999987E6 

1E6 

1000 

.100 

4.0E-09 

99. 

1.00000049999987E6 

1E6 

1000 

.010 

1.2E-10 

P0  =  .999 

=  1.0 

Case  No. 

R 

h 

k 

<Ty 

RERR 

100. 

3.33463215647498E0 

0.00 

0.00 

.500 

2.3E-11 

101. 

3.29205427406123E0 

0.00 

0.00 

.100 

8.0E-11 

102. 

3.29052673150727EO 

0.00 

0.00 

.010 

6.7E-16 
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P0  =  .999  Wx  =  1.0 


Case  No. 

R 

h 

k 

<Ty 

RERR 

103. 

3.53577235751072E0 

0.00 

1.00 

.500 

8.9E-16 

104. 

3.44218941775034E0 

0.00 

1.00 

.100 

8.0E-11 

105. 

3.43912287816784E0 

0.00 

1.00 

.010 

1.1E-07 

106. 

4. 1243 1141 2066 1 6E0 

1.00 

0.00 

.500 

7.4E-11 

107, 

4.09151266085973E0 

1.00 

0.00 

.100 

8.0E-11 

108. 

4.09028540157306E0 

1.00 

0.00 

.010 

1.1E-16 

109. 

6.93719921979876E0 

1.00 

5.00 

.500 

1.2E-15 

110. 

6.47656370770719E0 

1.00 

5.00 

.100 

8.0E-11 

111. 

6.46007535448830E0 

1.00 

5.00 

.010 

8.0E-11 

112. 

8.17462362868007E0 

5.00 

1.00 

.500 

3.9E-11 

113. 

8. 15266008069658E0 

5.00 

1.00 

.100 

8.0E-11 

114. 

8.15180095237486E0 

5.00 

1.00 

.010 

2.9E-08 

115. 

1.32590217041336E1 

5.00 

10.0 

.500 

2.2E-14 

116. 

1.28786495588648E1 

5.00 

10.0 

.100 

8.0E-11 

117.; 

1.28629674529177E1 

5.00 

10.0 

.010 

8.0E-11 

118. 

1.40783138161510E1 

10.0 

5.00 

.500 

7.7E-11 

119. 

1.40151628274090E1 

10.0 

5.00 

.100 

8.0E-11 

120. 

1.40126436416327E1 

10.0 

5.00 

.010 

9.1E-08 

121, 

5. 1 1535358789872E2 

100. 

500. 

.500 

6.9E-11 

122. 

5.10586815818419E2 

100. 

500. 

.100 

6.3E-11 

123. 

5.10517724887103E2 

100. 

500. 

.010 

8.3E-11 

124. 

5. 12947737652262E2 

500. 

100. 

.500 

8.0E-11 

125. 

5. 12933140365073E2 

500. 

100. 

.100 

8.0E-11 

126. 

5. 1293253 1471 364E2 

500. 

100. 

.010 

2.1E-08 

127, 

1.00000204511886E6 

1000. 

1E6 

.500 

7.8E-11 

128. 

1.00000080903890E6 

1000. 

1E6 

.100 

7.5E-11 

129. 

1 .00000053 1 05685E6 

1000. 

1E6 

.010 

2.1E-10 

130. 

1.00000359023064E6 

1E6 

1000. 

.500 

1.7E-09 

131. 

1.00000359023064E6 

1E6 

1000. 

.100 

6.9E-11 

132. 

1.00000359023064E6 

1E6 

1000. 

.010 

2.4E-13 
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ALGORITHM  AND  FORTRAN  LISTING  FOR  AERF 

In  this  appendix  the  algorithm  used  in  PKILL  for  the  computation  of  the  function  aerf  is 
given  where,  for  all  real  x  and  h, 

f(x, h)  =  aerf(x,  h)  =  erf(x  +  h)  —  erf(x  -  h)  A(l) 

or 

f(x,h)  =  aerf(x,  h)  =  erfc(x-h)  -erfc(x  +  h)  A(2) 

erf  y = i }  q  E(z) dz’  E(z) = exp(  - z2) 

o  f00 

erfc  y  =  1  -erf  y  =  j  E(z)  dz  . 

(As  an  aside,  note  that  i  ^  E(z)  dz  =  aerf  (x,  h)  where  x  =  — g—  and  h  =  --y-*,  ) 

It  is  sufficient  hereafter  to  consider  x  >  0  and  h  >  0,  since 

f(-x,h)  =  f(x,h),  f(x, -h)  =  -f(x,h)  and  f(0,  h)  =  2  erf  h,  f(x,0)  =  0. 

Also,  if  (x  -  h)  >  9.0,  then  f(x,h)  can  be  set  to  zero  if 

(x  -  h)2  +  2.76959  >  -  DXPARG(l).  A(3) 

Here  DXPARG(l)  is  the  smallest  negative  number  such  that 

exp[DXPARG(l)]  >  0. 

For  the  IBM  PC  double  precision  arithmetic  DXPARG(l)  ~  —708.389.  Also,  A(3)  follows  from 
f(x’h)<^j^E(z)  ^fh  dz  <  Jf  fx  -1)  <  exP  [  ~  (x  -  h)2  —  In  (9  )]  , 


For  sufficiently  small  h  and/or  x  some  approximations  are  used  that  improve  efficiency.. 


Consider 


and  let  u  =  (t  -  x)/h  to  obtain 


erf  (x  +  h)  =  erf  (x)  -f  j  X+^  E(t)  dt, 


Similarly, 


erf(x  +  h)  =  erf(x)  -f  E(x  +  hu)  du. 
erf(x-h)  =  erf  (x)  - J  *  E(x-hu)  du, 

f(x,h)  =  wJo  tE(x_hu)  +  E(x  +  hu)]  du> 

f  (x,  h)  =  E(x)  J  *  E(h  u)  cosh(2  x  h  u)  du  . 


A-l 
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If  h  and  x  are  both  sufficiently  small,  then  the  exponentials  in  A(5)  can  be  replaced  by  the  first 
few  terms  of  their  Maclaurin  expansions,  i.e, 

E(x  -  h)  +  E(x  +  h)  ~  1  -  (x  -  h)2  +  i(x  -  h)4  +  1  -  (x  +  h)2  +  £(x  +  h)4. 

Integration  yields 

f(x,h)  ~  §[(1  -x2  -Ih2)  +  I(x4  +  2x2h2  -f  |h4)]  . 


Thus 

provided 


f(x.h)^^(l-x2-ih2)  , 

|(T4  +  2T4 +iT4)  =  1.60  T4  <  EPS  , 


A(7) 


where  T  =  max(x,h),  and  EPS  is  given  by  (11). 

Another  approximation  that  is  used  can  be  obtained  from  (A6)  by  using  the  fact  that  if 
(2hx)2<2EPS,  then  cosh  2hxu  ~  1.0  , 


and 

f(x, h)  =  2  £(x)  erf  h  if  hx  <  ^EPS/2  .  A(8) 

In  addition,  if 

h2  <  3  EPS,  then  erf  h  ~  h, 

and 

f(x,h)~||E(x),  if  hx<4EP§72  and  h<>JTEPS.  A(9) 

The  final  approximation  used  is  given  by 

f(x,h)  ~erfc(x-h),  if  x  +  h>5.736,  withh>x,  or  if  exp(  — 4hx)  <  EPS.  A(10) 

In  the  first  case  erfc(x-h)  >  1,  and  erfc(x  +  h)  <  4.99  10  -16.  In  the  second  case,  we  have  by  using 

(A2)  that  of00 

f(x,h)=^j  E(u-h)[l-exp(-4hu)]  du. 

Hence,  since 


4hu  >  4hx  >  -log(EPS),  x  <  u  <  oo,  A(ll) 

the  approximation  given  by  (A  10)  holds.  A  subprogram  DESPLN,  written  by  A.  Morris  and  contained 
in  NSWCLIB,  supplies  log(EPS),  where  EPS  is  given  in  (11) . 

In  case  none  of  the  approximations  above  are  applicable  for  evaluating  aerf,  then  two  other 
procedures  are  used  which  follow.  One  is  used  if  x  <  0.4  and  the  other  is  used  when  x  >  0.4  . 


If  x  <  0.4 

f(x,  h)  = 

1 

This  result  follows  directly  by  expanding  (2/V?)E(z)  into  its  Maclaurin  series  and  integrating  the  series 
from  x  —  h  to  x  +  h. 


g  (_  i)"  [(x  +  h)2n  +  1-(x -h)2n  +  1] 

i  =  0 


n! 


2n  + 1 


A(12) 


A-2 
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For  the  computation  of  f(x,  h)  in  PKILL  let 


Sn  — 

uu_v  =uS„-i  +  v"  \  n>1»  Sj  =  1. 

A(13) 

Then 

f(x,h) 

_.4h^  (-1)"  , 

(2n  +  1)  2n  +  1’ 

A(14) 

where,  from  (A12)  and  (A13), 

^2n  +  l 

=  (x  +  h)  S2n  +  (x  -  h)2” 

=  (x  +  h)[(x  +  h)S2„_1  +  (x-h)2n-1]  +  (x-h)2" 

=  (x  +  h)2  S2n_1 +  2x(x-h)2n_1  ,  n  >  1.  A(15) 

If  x  >  0.40 

f(x,  h)  =  ^  E(x/tf )  £e(x/>[2  )  ,  A(16) 

where  Hn  denotes  the  Hermite  polynomial  of  degree  n  [1,  pp.  781-802].  This  series  for  aerf  follows  by 
obtaining  the  Maclaurin  expansions  of  erf(x  +  h)  and  erf(x-h)  and  subtracting  the  resulting  series. 


The  expansions  for  erf(x  ±  h)  are  derived  by  using  the  following  results: 

^[erfz]=|=E(z)  A(  17.1) 

£B[E(z)]  =  (-l)"E(z)Hn(z)  A(17.2) 

[erfz]  =  ( -  1)"  $  E(z)  Hn(z)  A(17.3) 

Hn  +  1(z)  =  2zHn(z)-2nHn_1(z),  H0  =  l,  fl1  =  2z,  A(17.4) 


and 

erf  (x  +  h)  =  erf  (x)  +  ^  E(x)  f)  •  A(18) 


For  the  computation  of  (A  16)  in  PKILL,  we  rewrite  (A16)  in 

the  form 

where 

,  00 
f(x,h)=|E(#)£ 

1  n  =  0 

J2n 

2n  +  1  ’ 

A(19) 

In  addition, 

n  >  0. 

A(20) 

^2n  =  [2hx  J2n-1  —  2h2  J2n  —  2 1/^ n  » 

n  >  1, 

J0  =  E(x/^)  , 

^2n  —  1  =  [2hx  J2n_2  — 2h2  ^2n  —  —  » 

n  >  2, 

J1  =  2hxE(x/^)  . 

A(21) 

A-3 
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Since 

H2n(x)  <  exp(x2/2)  22n  + 1  n!  ,  [1,  p.  787]  A(22) 


use  of  Stirling’s  formula  for  the  factorial  gives  an  estimate  for  J2„  ,  i.e., 


Now  since  A(19)  is  only  used  when  4hx  <  -  log(EPS)  and  x  >  3h  (see  A(ll)  and  page  A-5),  it  follows 
even  for  a  machine  with  a  19  digit  mantissa  that  h2  <  43.75/12  ~  3.6458.  Hence,  using  A(23),  no  more 
than  28  terms  of  A(19)  would  be  needed  to  maintain  a  relative  error  less  than  10“ 12  in  f(x,h)  . 
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DOUBLE  PRECISION  FUNCTION  AERF  (X,H) 

C - - - 

C  COMPUTATION  OF  ERF(X  +  H)  -  ERF(X  -  H) 

C- - 

DOUBLE  PRECISION  ERF  ,ERFC  ,DEPSLN  .DPMPAR  ,DXPARG 
DOUBLE  PRECISION  AH  ,AX  ,A2  ,C  ,C1  ,DN2 

DOUBLE  PRECISION  E  ,EPS  ,H  ,HF  ,HG  ,H2 

DOUBLE  PRECISION  H3  ,S  ,ST  ,T  ,U  ,V  ,X 
DOUBLE  PRECISION  XMH  ,XPH  ,X2  ,Z 

C- - 

C  C  =  2/DSQRT(PI),  Cl  =  In  (9*DSQRT(PI)) 

C- - 

DATA  C  /1.12837916709551257D0/,  C1/2.76959D0/ 

C - 

C  EPS  IS  A  MACHINE  DEPENDENT  CONSTANT.  EPS  IS  THE  SMALLEST 
C  DOUBLE  PRECISION  FLOATING  POINT  NUMBER  FOR  WHICH 
C  1.0  +  EPS  .GT.  1.0. 

C- - 

EPS  =  DPMPAR(l) 

C - - 

AERF  =  0.D0 

IF  (H  .EQ.  0.D0)  RETURN 
C 

AH  =  DABS(H) 

AX  =r  DABS(X) 

XPH  =  AX  +  AH 
XMH  =  AX  -  AH 
C 

T  =  DMAX1(AH,  AX) 

T  =  T*T 

IF  (  1.6*T*T  .LT.  EPS)  GOTO  140 
IF  ((AH*AX)**2  .LE.  EPS/2)  GOTO  150 
IF  (AX  .LE.  AH)  GOTO  100 
C 

IF  (XMH  .LT.  9 DO)  GOTO  5 

IF  (XMH*XMH  +  Cl  .GT.  -DXPARG(l) )  RETURN 
5  IF  (4*AH*AX  .GT.  -DEPSLN(O))  GOTO  120 
C 

IF  (AX  .GT.  3D0*AH)  GOTO  10 
IF  (XPH  .LT.  1.D0)  GOTO  110 
GOTO  130 
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FOR  (AX  LESS  THAN  OR  EQUAL  TO  .40) 


10  E  =  DMAX1(1.D-15,  EPS) 

IF  (AX  .GT.  0.40D0)  GOTO  30 

H2  =  XPH+XPH 

A2  =  XMH+XMH 

X2  =  AX  +  AX 

ST  =  l.DO 

HF  =  XMH 

N  =  0 

N1  =  1 

DN2  =  l.DO 

S  =  O.DO 

20  N  =  N  +  1 
N1  =  N1  +  2 
DN2  =  -DN2/N 
ST  =  H2*ST  +  X2*HF 
HF  =  A2*HF 
T  =  ST+DN2/N1 
S  =  S  +  T 

IF  (DABS(T)  .GT.  E*DABS(S))  GO  TO  20 
S  =  0.5D0  +  (S  +  0.5D0) 

AERF  =  2*C*AH*S 
GO  TO  45 


FOR  (AX  GREATER  THAN  .40) 


N  =  1 
J  =  0 
H2  =  O.DO 

Z  =  DEXP(-AX*AX/2) 
U  =  2*AH*C*Z 
H3  =  Z 
V  =  2*H*H 
HF  =  2*AX*AH 


S  =  O.DO 

35  H2  =  (HF*H3  -  V*H2)/N 
N  =  N  +  1 

H3  =  (HF*H2  -  V*H3)/N 
N  =  N  +  1 
HG  =  H3/N 
S  =  S  +  HG 

IF  (DABS(HG)  .GT,  E*DABS(S))  GO  TO  35 
IF  (J  .NE.  0)  GO  TO  40 
J  =  1 
GO  TO  35 

40  AERF  =  U*(S  +  Z) 

45  IF  (H  .LT.  O.DO)  AERF  = -AERF 
RETURN 
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C- - 

C  SPECIAL  CASES 
C - 

100  IF  (XPH  .LT.  5.736D0)  GOTO  110 
IF  (XMH  .GT.  -5.6D0)  GOTO  120 
AERF  =  SIGN(2.D0,  H) 

RETURN 

C 

110  AERF  =  ERF(XPH)  -  ERF(XMH) 

IF  (H  .LT.  0.D0)  AERF  = -AERF 
RETURN 
C 

120  AERF  =  ERFC(XMH) 

IF  (H  .LT.  0.D0)  AERF  = -AERF 
RETURN 
C 

130  AERF  =  ERFC(XMH)  -  ERFC(XPH) 

IF  (H  .LT.  0.D0)  AERF  = -AERF 
RETURN 
C 

140  AERF  =  2*C*H*(.5D0  +  (.5D0  -  (X*X  +  H*H/3))) 
RETURN 
C 

C  THE  VALUE  IS  2*DEXP(-X*X)*ERF(H) 

C 

150  T  =  2.D0 
X2  =  X*X 

IF  (X2  .GE.  EPS)  T  =  2*DEX?(-X*X) 

IF  (H*H  .GE.  3*EPS)  GOTO  160 
AERF  ss  C*H*T 
RETURN 

160  AERF  =  T*ERF(H) 

RETURN 

END 
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APPENDIX  B 

FORTRAN  LISTING  FOR  PKILL 
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SUBROUTINE  PRILL  (R0,SX,SY,H,R,P) 

C  PRILL  GIVES  THE  ELLIPTICAL  COVERAGE  FUNCTION,  WHERE  P  DE- 
C  NOTES  THE  PROBABILITY  OF  A  SHOT,  NORMALLY  DISTRIBUTED  WITH 
C  MEAN  (0, 0)  AND  STANDARD  DEVIATIONS  SX,  SY  IN  THE  X,Y 
C  DIRECTIONS,  RESPECTIVELY,  FALLING  IN  A  CIRCLE  IN  THE  XY- 
C  PLANE  OF  RADIUS  R0  AND  CENTERED  AT  (H,  R).  THE  INPUT  IS  R0, 

C  SX,  SY,  H,  R.  THE  OUTPUT  IS  P. 

C  PRILL  GIVES  AT  LEAST  6-SIGNIFICANT  DIGIT  ACCURACY  FOR  P  PRO- 
C  VIDED  MAX(H/S,R/S,MAX(SX,SY)/S)  .LE.  10/DSQRT(DPMPAR(1)),  WHERE 
C  S  =  MIN(SX,SY).  IF  P  .LT.  l.E-20  ACCURACY  MAY  BE  LESS  THAN  6  SD. 

C  IF  SX  OR  SY  .LE.  0,  THEN  P  SET  TO  ( -  1E10). 

C 

C  REFERENCES... 

C  NWL  REPORT  N0.1710,  AUG.1960.  MATH  OF  COMP.  OCT.  1961, 

C  PP.  375,382.  NSWC  REPORT  NO.83-13,  NOVEMBER,  1982. 

C- - 

DOUBLE  PRECISION  R,  R8,  V(16),  W(16) 


V(*),  W(*)-  GAUSSIAN  ABSCISSAS  AND  WEIGHTS  OF  ORDER  32,  ON  (-1,1). 


EXTERNAL  AERF  , 

DPMPAR  ,DXPARG  ,DEPSLN  ,ERF  ,ERFC 

,GRATIO 

DOUBLE  PRECISION 

A 

AERF  ,A1 

,A2  ,A3  ,A4 

DOUBLE  PRECISION 

C 

,C8 

,D  ,DEPSLN  ,DM  .DPMPAR 

,DXPARG 

DOUBLE  PRECISION 

D1 

,D2 

,EPS 

,ERF  ,ERFC  ,E3 

DOUBLE  PRECISION 

E4 

,F 

,G2 

,G3  ,H  ,H1 

DOUBLE  PRECISION 

H2 

,H3 

,H5 

,H6  ,H8  ,P 

DOUBLE  PRECISION 

Q 

,Qi 

,R 

,RPI  ,RT2  ,R0 

DOUBLE  PRECISION 

R8 

,SA 

,SEPS1  ,SQEPS  ,SX  ,SY 

DOUBLE  PRECISION 

SO 

,S1 

,S2 

,S9  ,T  ,T1 

DOUBLE  PRECISION 

T2 

,T3 

,T4 

,T5  ,T9  ,U 

DOUBLE  PRECISION 

XI 

,X5 

,YY 

,Z  ,ZA  ,Z1 

DOUBLE  PRECISION 

Z2 

,Z3 

,Z8 

,Z9 

DATA  V(l)  /.4830766568773832D-01/,  V(2)  /.1444719615827965D+00/, 

*  V(3)  /.2392873622521371D+00/,  V(4)  /.3318686022821276D+00/, 

*  V(5)  /.4213512761306353D+00/,  V(6)  /.5068999089322294D+00/, 

*  V(7)  /.5877157572407623D+00/,  V(8)  /.6630442669302152D+00/, 

*  V(9)  /.7321821187402897D+00/,  V(10) /.7944837959679424D+00/, 

*  V(ll)  /.8493676137325700D+00/,  V(12)  /.8963211557660521D+00/, 

*  V(13)  /.9349060759377397D+00/,  V(14)  /.9647622555875064D+00/, 

*  V(15)  /.98561 151 15452683D+00/,  V(16)  / .99726386 1 84948 1 6D-I-00/ 
DATA  W(l)  /.9654008851472780D-01/,  W(2)  /.9563872007927486D-01/, 

*  W(3)  /.9384439908080457D-01/,  W(4)  /.9117387869576388D-01/, 

*  W(5)  /.876520930044038 ID-01  /,  W(6)  /.8331192422694676D-01/, 

*  W(7)  /.7819389578707031D-01/,  W(8)  /.7234579410884851D-01/, 

*  W(9)  /.6582222277636185D-01/,  W(10)  /.5868409347853555D-01/, 

*  W(ll)  /.5099805926237618D-01  /,  W(12)  /.4283589802222668D-01/, 

*  W(13)  /.3427386291302143D-01/,  W(14)  /.2539206530926206D-01/, 

*  W(15)  /.  1627439473090567D-01/,  W(16)  /.7018610009470097D-02/ 
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C  =  (2**.25)*GAMMA(3/4)/PI 
RPI  =  1/DSQRT(PI) 

RT2  =  DSQRT(2) 

DATA  A3  /20.0D0/ 

DATA  C  /.46386480428950042D0/ 

DATA  RT2  /1.414213562373095D0/ 

DATA  RPI  /.5641895835477563D0/ 

P  =  0.0D0 

IF  (SX  .GT.  0.0D0  .AND.  SY  .GT.  0.0D0)  GOTO  5 
P  =  -1.D10 
RETURN 
5  J  =  0 

EPS  =  DPMPAR(l) 

SQEPS  =  DSQRT(EPS) 

A  =  6.5D0 


TEST  NO.l  (REPORT  83-13). 


Z3  =  DMIN1(DSQRT(DSQRT(DPMPAR(2))),  l.D-30)*SX*SY 
IF  (R0+R0  .LE.  Z3)  RETURN 


H2  =  H*H  +  K*K 
H8  =  DABS(H) 

K8  =  DABS(K) 

DM  =  DMAX1(SX,SY) 


TEST  N0.2  (REPORT  83-13). 


IF  ((R0  -  H8  +  A3*SX)  .LE.  0.0D0)  RETURN 
IF  ((R0  -  K8  +  A3*SY)  .LE.  0.0D0)  RETURN 

C- - 

C  TEST  N0.3  (REFORT  83-13). 

C- - 

C8  =  -DEPSLN(O) 

C- - DEXP(- 10.3026)  =  3.35E(-5) - 

C- - C8  =  -DLOG(EPS) - 

A4  =  DMAX1(C8  -  1.74249D1,  10.3026D0) 

A4  =  DSQRT(A4  +  A4) 

H3  =  (R0  -  H)*(R0  +  H) 

H5  =  (R0  -  K)*(R0  +  K) 

T  =  R0  -  A4+DM 

IF  (T  .LT.  0.0D0)  GOTO  25 

IF  (T*T  .LT.  H2)  GOTO  25 

IF  (R0  .LT.  1/SQEPS)  GOTO  20 

T  =  H3  -2.0D0*R0*A4*DM  +  (A4*DM)**2  -  K*K 
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IF  (DABS(H3)  .LE.  DABS(H5))  GOTO  10 
T  =  H5  -2.0D0*R0*A4*DM  +  (A4*DM)**2  -  H*H 
10  IF  (T  .LT.  0.0D0)  GOTO  25 
20  P  =  1.0D0 
RETURN 

C  TEST  N0.4  (REPORT  83-13). 

25  SO  =  DSQRT(H2) 

G2  =  0.D0 

IF  (SO  .LE.  R0)  GOTO  30 
D  =  ((SO  -  R0)/DM)**2 

IF  (RO*RO*DEXP(-O.5D0*D)  .LE.  Z3)  RETURN 

C  SX  -  SY  .LT.  EPS 

C- - 

30  IF  (DABS(SX  -  SY)  .GT.  20.0D0*DM*EPS)  GOTO  35 
118  =  SO 
K8  =  0.0D0 

IF  ((R0  -  H8  +  A3*DM)  .LE.  0.0D0)  RETURN 
IF  (R0  .LT.  DM/SQEPS)  GOTO.  35 
J  =  1 

G2  =  (K*K  -  H3)/((R0  +  H8)*DM) 

IF  (DABS(HS)  .LE.  DABS(H5))  GOTO  35 
G2  =  (H*H  -  H5)/((R0  +  H8)*DM) 

C  NORMALIZE  DMAXl(SX.SY)  =  1. 

35  R  =  R0/DM 

51  =  SX/DM 

52  =  SY/DM 
H8  =  H8/DM 
K8  =  K8/DM 

H2  =  H2/(DM*DM) 

C- - 

C  SI  =  1  .GE.  S2 

C- - 

IF  (SI  .GE.  S2)  GOTO  40 
HI  =  SI 

51  =  S2 

52  =  HI 
HI  =  H8 
H8  =  K8 
K8  -  HI 

C- - 

C  LIMITING  RESULTS  FOR  DMIN1(S1,S2)  =  0 
C- - 
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40  SEPS1  =  DMIN1(6.71D0*SQEPS,  l.D-5) 


R  =  K,  S2  SMALL 


IF  (K8  .NE.  R)  GOTO  45 

YY  =  .166484D0*(R*(H8*H8  +  1.0D0)  +  1.0D0/R)*S2 
IF  (DABS(YY)  .GT.  SEPS1)  GOTO  85 
HI  =  H8/RT2 

P  •■=  C*DEXP(-H1*H1)*DSQRT(K8*S2) 

RETURN 


R  .GT.  K,  S2  SMALL 


45  IF  (R  .LT.  K8)  GOTO  85 
IF  (K8  .NE.  0.0D0)  GOTO  75 
HI  =  S2*S2/(4.0D0*RT2*R) 

G2  =  G2/RT2 

IF  (J  .EQ.  0)  G2  =  DABS(H8  -  R)/RT2 
IF  (DABS(G2)  .LT.  4.0D0)  GOTO  55 
50  IF  (H1*DABS(G2)  .GT.  SEPS1)  GOTO  85 
P  =  .5D0*AERF(H8/RT2,R/RT2) 

RETURN 

55  IF  (J  .NE.  0)  GOTO  60 

P  =  .5D0*AERF(H8/RT2,R/RT2) 

GOTO  70 

60  IF  (H8  +  R  .LT.  RT2)  GOTO  65 

P  =  .5D0*(ERFC(G2)  -  ERFC((H8  +  R)/RT2)) 

GOTO  70 

65  P  =  ,5D0*(ERF((H8  +  R)/RT2)  -  ERF(G2)) 

70  IF  (H1*DEXP(-G2*G2)  .GT.  P*SEPS1)  GOTO  85 
RETURN 

75  Z  =  (R  -  K8)*(R  +  K8) 

G2  =  DSQRT(Z) 

HI  =  DABS(H8  -  G2) 

J  =  0 

IF  (HI  .GT.  5.D0)  GOTO  80 
J  =  1 

Z8  =  AERF(H8/RT2,G2/RT2) 

IF  (Z8  .EQ.  0.0D0)  GOTO  85 
HI  =  DEXP(-0.5DQ*(H8  -  G2)**2)/Z8 

80  HI  =  1.0D0/(8.0D0*RT2*Z)*(K8*K8*DABS(H8  -  G2)  +  R^R/G2)*S2*S2*H1 
IF  (DABS(Hl)  .GT.  SEPS1)  GOTO  85 
U  =  .5D0 

IF  (J  .EQ.  0)  Z8  =  AERF(H8/RT2,G2/RT2) 

S2  =  RT2*S2 

IF  (K8  -  R  .GT.  -13.D0*S2)  U  =  .25D0*AERF(K8/S2,R/S2) 

P  =  U*Z8 
RETURN 
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C - 

C  FIND  THE  VALUE  OF  A 


85  J8  =  0 
SO  =  Sl*Sl 
S9  =  S2*S2 
G3  =  H8+H8 
G2  =  K8*K8 
Z8  =  SO  +  S9 
Z  =  S0*S0  +  S9*S9 
H3  =  G3+S0  +  G2+S9 
T1  =  2.0D0*(Z  +  2*H3) 

YY  =  R*R*(H2  +  Z8)/T1 
T1  =  (H2  +  Z8)*(H2  +  Z8)/T1 
CALL  GRATIO(T1,YY,Z,Z8,0) 

Z2  =  R/(RT2*S2) 

R8  =  R/(RT2*S1) 

H2  =  H8/(RT2*S1) 

H3  =  K8/(RT2*S2) 

SO  =  0.D0 
S9  =  0.D0 

IF  (Z  .LT.  l.D-13  .AND.  S2  .GT.  5.D-10)  GOTO  90 
IF  (H2  .GT.  50.0D0  .OR.  H3  .GT.  50.0D0)  SO  =  1.5D0 
90  U  =  AERF(H3,Z2) 

YY  =  .25D0*U*AERF(H2,R8) 

IF  (YY  .GT.  0.ID0)  S9  =  .5D0 
IF  (YY  ,GE.  5.D-15)  GOTO  95 
Z  =  YY 
GOTO  !00 

95  Z  =  DMIN1(YY,Z) 

100  IF  (Z  .GE.  5.D-1)  GOTO  105 
A  =  A  +  .5D0 

IF  (Z  .GE.  5.D-4)  GOTO  105 

A  =  A  +  .5D0  +  SO 

IF  (Z  .GE.  l.D-6)  GOTO  105 

A  =  A  +  .5D0  +  S9 

IF  (Z  .GE.  5.D-9)  GOTO  105 

A  =  A  +  .5D0 

IF  (Z  .GE.  5.D-10)  GOTO  105 
A  =  A  +  .25D0 

IF  (Z  .GE.  5.D-11)  GOTO  105 
A  =  A  +  .25D0 

IF  (Z  .GE.  5.D-I2)  GOTO  105 
A  =  A  +  .25D0 

IF  (Z  .GE.  5.D-13)  GOTO  105 

A  =  A  +  .5D0  +  S9/2 

IF  (Z  .GT.  5.D-15)  GOTO  105 
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A  =  A  +  .50D0 

IF  (Z  .GT.  l.D-18)  GOTO  105 
A  =  A  +  .25D0 

IF  (Z  .GT.  l.D-20)  GOTO  105 
A  =  A  +  .25D0 

IF  (Z  .GT.  l.D-25)  GOTO  105 
A  =  A  +  1.D0 

IF  (Z  .GT.  l.D-30)  GOTO  105 
A  =  A  +  2.D0 

105  IF  (S2  .GE.  5D-2  .OR.  H3  .LE.  Z2)  GOTO  107 
T9  =  R8*U*DEXP(-H2*H2) 

IF  (T9  .LT.  5D-2+Z)  A  =  A  +  .5D0 

C - 

C  START  INTEGRATION  PROCEDURE 

107  SO  =  SI 
S9  =  S2 
Z8  =  S2 
G2  =  K8 
G3  =  H8 
Z9  =  S1 

C - 

C  DETERMINE  INTERVAL  OF  INTEGRATION,  (A,B). 
C  E3  =  (B-A)/2,  D1  =  (B+A)/2. 

110  Z  =  G2  +  A*Z8 
H3  =  G2  -  A*Z8 
H5  =  0.0D0 
T3  =  -1.0D0 
A1  =  (G3  -  A*Z9)/R 

IF  (DABS(A1  -  0.5D0)  .GE.  0.5D0)  GOTO  115 
A2  =  ((1.0D0  -  G3/R)  +  A*Z9/R)*(1.0D0  +  Al) 

SA  =  DSQRT(A2) 

IF  (SA  .LT.  H3/R)  GOTO  115 
T3  =  A1/DSQRT(1.0D0  +  SA) 

C - 

C  T9  =  1.0  -  D1 

C- - 

115  IF  (H3  .GT.  0.0D0)  GOTO  135 

C - 

C  H3  .LE.  0.0 

IF  (T3  .LT.  0.0D0)  GOTO  120 
D2  =  0.5D0*(1.0D0  +  T3) 

E4  =  0.25D0*SA/D2 
T5  =  E4 

120  IF  (Z  .LT.  R)  GOTO  125 


B-6 


QOQ 


NAVSWC  TR  90-513 


E3  =  .5D0 
D1  =  E3 
T9  =  D1 
GO  TO  130 

125  D1  =  .5D0*(1.0D0  +  DSQRT(1.0D0  -  Z/R)) 
E3  =  .25D0*Z/(R*D1) 

T9  =  E3 

130  IF  (T3  .LE.  D1  -  E3)  GOTO  165 
GOTO  160 

H3  .GT.  0.0 


135  H5  =  l.ODO 
Q  =  H3/R 

IF  (Q  .LE.  l.ODO)  GOTO  140 
E3  =  .5D0 
D1  =  E3 
T9  =  D1 
GOTO  165 

140  IF  (T3  .LT.  0.0D0)  GOTO  145 

E4  =  DSQRT((1.0DO  -  G2/R)  +  A*Z8/R) 

D2  =  0.5D0*(E4  +  T3) 

T5  =  0.5D0*(Q/(1.D0  +  E4)  +  SA/(1.D0  +  T3)) 

E4  =  0.25D0*(SA  -  Q)/D2 
145  IF  (Z  .LT.  R)  GOTO  150 

E3  =  0.5D0*DSQRT((1.0D0  -  G2/R)  +  A*Z8/R) 

D1  =  E3 

T9  =  0-25D0*(3.0D0  +  H3/R)/(1.0D0  +  E3) 

GOTO  155 

150  T9  =  DSQRT((1.0D0  -  G2/R)  +  A+Z8/R) 

T2  =  DSQRT(1.0D0  -  Z/R) 

D1  =  0.5D0*(T9  +  T2) 

E3  =  A*Z8/(R*2.0DO*D1) 

T9  =  O.5D0*((H3/R)/(l.ODO  +  T9)  +  (Z/R)/(1.0D0  +  T2)) 
155  IF  (T3  .LE.  D1  -  E3)  GOTO  165 
160  E3  =  E4 
D1  =  D2 
T9  =  T5 

165  IF  (J8  .NE.  0)  GOTO  170 
J8  =  1 
F  =  E3 
T  =  D1 
T1  =  T9 
H6  =  H5 
Z8  =  SI 
G2  =  H8 
G3  =  K8 
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Z9  =  S2 
GOTO  110 


DETERMINE  IN  WHICH  ORDER  THE  X  AND  Y  INTEGRATIONS  ARE  CARRIED  OUT. 


170  IF  (S2  .GT.  2D-2  .AND.  H8  +  K8  .LT,  2D2)  GOTO  172 
IF  (DABS(E3  -  F)  .GT.  .4D0*F)  GOTO  172 
IF  (D1  -  T)  200,  180,  195 
172  IF  (E3  .LT.  2.D4*EPS)  GOTO  195 
IF  (F  .LT.  2.D4*EPS)  GOTO  200 
IF  (DMAX1(H8/S1,K8/S2)  .LT.  2.0D0)  GOTO  175 
IF  (S2  .LT.  l.D-5)  GOTO  175 
IF  (S2  .GE.  5.D-4)  GOTO  180 
IF  (SI  .NE.  S2)  GOTO  195 

175  IF  (DMIN1(E3,  F)  .LT.  2.5D-2+SQEPS)  GOTO  185 
180  IF  (E3  -  F)  200,  190,  195 
185  IF  (E3  -  F)  195,  190,  200 
190  IF  (SO  .GE.  S9)  GOTO  200 
195  E3  =  F 
D1  =  T 
T9  -  T1 
S9  =  SI 
SO  =  S2 
Z8  =  H8 
H8  =  K8 
K8  =  Z8 
H5  =  H6 

200  Z2  =  R/(RT2*S9) 

R8  =  R/(RT2*S0) 

H2  =  H8/(RT2*S0) 

H3  =  K8/(RT2*S9) 

N1  =  16 
203  IZ  =  0 
IZ3  =  0 
IY  =  0 
P  =  0.D0 

T  =  H2  -  R8  +  R8*(D1  -  E3*V(16))**2 
IF  (T  .GT.  0.D0  .AND.  T*T  .GT.  -DXPARG(l))  RETURN 
Q1  =  RPI*E3*R8 
G3  =  0.0D0 
NT  =  2*N1  +  1 
Z  =  (.5D0  +  Dl)  +.5D0 
DO  260  II  =  1,  NT 
205  I  =  II  -  (N1  +  1) 

IF  (I  .EQ.  0)  GOTO  260 
J  =  IABS(I) 

Q  =  E3*ISIGN(1,I)*V(J) 
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T  =  Q  +  D1 
H6  =  Z  +  Q 
Q  =  T9  -  Q 
F  =  H6*Q 
T1  =  R8*F 

T2  =  (H2  -  T1)*(H2  -  Tl) 

IF  (H2  -  R8  .GE.  0.D0  .OR.  T  .LT.  EPS)  T2  = 

*  ((H2  -  R8)  +  R8*T*T)**2 

T4  =  DEXP(-T2) 

IF  (H2  .NE.  0.0D0)  GOTO  210 
T4  =  T4  +  T4 
GO  TO  215 

210  IF  (H5  .NE.  0.0D0)  GOTO  215 
T2  =  4.0D0*H2*T1 
IF  (T2  .GT.  C8)  GO  TO  215 
T4  =  T4*(1.0D0  +  DEXP(-T2)) 

215  IF  (IZ  .NE.  0)  GOTO  255 
G2  =  DSQRT(1.0D0  +  F) 

Z1  =  Z2*T*G2 
XI  =  H3  -  Z1 

IF  (XI  .GT.  -A)  GOTO  225 
220  IZ  =  1 

X5  =  2.D0 
GO  TO  255 

225  IF  (DABS(Xl)  .GT.  l.D-2*Zl)  GOTO  230 
IY  =  1 

XI  =  ((K8  -  R)  +  R*F*F/(1.D0  +  T*G2))/(RT2*S9) 

230  IF  (IZ3  .NE.  0)  GOTO  250 

SA  =  H3  +  Z1 

IF  (SA  .GT.  A3)  GOTO  245 
IF  (IY  .EQ.  0)  GOTO  240 
IF  (XI  .GT.  RT2)  GOTO  235 
X5  =  ERF(SA)  -  ERF(Xl) 

GOTO  255 

235  X5  =  ERFC(Xl)  -  ERFC(SA) 

GOTO  255 

240  X5  =  AERF(H3,Z1) 

GOTO  255 
245  IZ3  =  1 

250  X5  =  ERFC(Xl) 

255  G3  =  G3  +  X5*T4*T*W(J) 

260  CONTINUE 
P  =  Q1*G3 

IF  (P  .GT.  YY)  P  =  YY 

IF  (P  .GT.  (1.0D0  -  DMIN1(1.D6*EPS,  l.D-5)))  P  =  1.0D0 

IF  (P  .LT.  0.0D0)  P  =  0.0D0 

RETURN 

END 
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